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The  Eshelby-type  problem  of  an  arbitrary-shape  polyhedral  inclusion  embedded  in  an 
infinite  homogeneous  isotropic  elastic  material  is  analytically  solved  using  a  simplified 
strain  gradient  elasticity  theory  (SSGET)  that  contains  a  material  length  scale  para¬ 
meter.  The  Eshelby  tensor  for  a  polyhedral  inclusion  of  arbitrary  shape  is  obtained  in  a 
general  analytical  form  in  terms  of  three  potential  functions,  two  of  which  are  the  same 
as  the  ones  involved  in  the  counterpart  Eshelby  tensor  based  on  classical  elasticity. 
These  potential  functions,  as  volume  integrals  over  the  polyhedral  inclusion,  are 
evaluated  by  dividing  the  polyhedral  inclusion  domain  into  tetrahedral  duplexes,  with 
each  duplex  and  the  associated  local  coordinate  system  constructed  using  a  procedure 
similar  to  that  employed  by  Rodin  (1996.  J.  Mech.  Phys.  Solids  44,  1977-1995).  Each  of 
the  three  volume  integrals  is  first  transformed  to  a  surface  integral  by  applying  the 
divergence  theorem,  which  is  then  transformed  to  a  contour  (line)  integral  based  on 
Stokes’  theorem  and  using  an  inverse  approach  different  from  those  adopted  in  the 
existing  studies  based  on  classical  elasticity.  The  newly  derived  SSGET-based  Eshelby 
tensor  is  separated  into  a  classical  part  and  a  gradient  part.  The  former  contains 
Poisson's  ratio  only,  while  the  latter  includes  the  material  length  scale  parameter 
additionally,  thereby  enabling  the  interpretation  of  the  inclusion  size  effect.  This  SSGET- 
based  Eshelby  tensor  reduces  to  that  based  on  classical  elasticity  when  the  strain 
gradient  effect  is  not  considered.  For  homogenization  applications,  the  volume  average 
of  the  new  Eshelby  tensor  over  the  polyhedral  inclusion  is  also  provided  in  a  general 
form.  To  illustrate  the  newly  obtained  Eshelby  tensor  and  its  volume  average,  three 
types  of  polyhedral  inclusions  -  cubic,  octahedral  and  tetrakaidecahedral  -  are 
quantitatively  studied  by  directly  using  the  general  formulas  derived.  The  numerical 
results  show  that  the  components  of  the  SSGET-based  Eshelby  tensor  for  each  of  the 
three  inclusion  shapes  vary  with  both  the  position  and  the  inclusion  size,  while  their 
counterparts  based  on  classical  elasticity  only  change  with  the  position.  It  is  found  that 
when  the  inclusion  is  small,  the  contribution  of  the  gradient  part  is  significantly  large 
and  should  not  be  neglected.  It  is  also  observed  that  the  components  of  the  averaged 
Eshelby  tensor  based  on  the  SSGET  change  with  the  inclusion  size:  the  smaller  the 
inclusion,  the  smaller  the  components.  When  the  inclusion  size  becomes  sufficiently 
large,  these  components  are  seen  to  approach  (from  below)  the  values  of  their  classical 
elasticity-based  counterparts,  which  are  constants  independent  of  the  inclusion  size. 
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1.  Introduction 

Eshelby’s  (1957,  1959)  solution  for  the  problem  of  an  infinite  homogeneous  isotropic  elastic  material  containing  an 
ellipsoidal  inclusion  prescribed  with  a  uniform  eigenstrain  is  a  milestone  in  micromechanics.  The  solution  for  the  dynamic 
Eshelby  ellipsoidal  inclusion  problem  was  obtained  by  Michelitsch  et  al.  (2003),  which  reduces  to  the  Eshelby  solution  in 
the  static  limiting  case.  Both  of  these  solutions  are  based  on  classical  elasticity.  Recently,  the  Eshelby  ellipsoidal  inclusion 
problem  was  solved  by  Gao  and  Ma  (2010b)  using  a  simplified  strain  gradient  elasticity  theory,  which  recovers  Eshelby’s 
(1957,  1959)  solution  when  the  strain  gradient  effect  is  not  considered. 

A  remarkable  property  of  Eshelby’s  (1957)  solution  is  that  the  Eshelby  tensor,  which  is  a  fourth-order  strain 
transformation  tensor  directly  linking  the  induced  strain  to  the  prescribed  uniform  eigenstrain,  is  constant  inside  the 
inclusion.  However,  this  property  is  true  only  for  ellipsoidal  inclusions  (and  when  classical  elasticity  is  used),  which  is 
known  as  the  Eshelby  conjecture  (e.g.,  Eshelby,  1961 ;  Rodin,  1996;  Markenscoff,  1998a, b;  Lubarda  and  Markenscoff,  1998; 
Liu,  2008;  Li  and  Wang,  2008;  Gao  and  Ma,  2010b;  Ammari  et  al.,  2010). 

For  non-ellipsoidal  polyhedral  inclusions,  Rodin  (1996)  provided  an  algorithmic  analytical  solution  and  showed  that 
Eshelby’s  tensor  cannot  be  constant  inside  a  polyhedral  inclusion,  thereby  proving  the  Eshelby  conjecture  in  the  case  of 
polyhedral  inclusions.  The  expressions  of  Eshelby’s  tensor  for  two-dimensional  (2-D)  polygonal  inclusions  were  included 
in  Rodin  (1996).  The  explicit  expressions  of  the  Eshelby  tensor  for  three-dimensional  (3-D)  polyhedral  inclusions  were 
later  derived  by  Nozaki  and  Taya  (2001),  where  an  exact  solution  for  the  stress  field  inside  and  outside  an  arbitrary-shape 
polyhedral  inclusion  was  obtained  and  numerical  results  for  five  regular  polyhedral  inclusion  shapes  and  three  other 
shapes  of  the  icosidodeca  family  were  presented.  Both  Rodin  (1996)  and  Nozaki  and  Taya  (2001 )  made  use  of  an  algorithm 
developed  by  Waldvogel  (1979)  for  evaluating  the  Newtonian  (harmonic)  potential  over  a  polyhedral  body.  A  more 
compact  form  of  the  Eshelby  tensor  than  that  presented  in  Nozaki  and  Taya  (2001 )  for  a  polyhedral  inclusion  in  an  infinite 
elastic  space  was  proposed  by  Kuvshinov  (2008)  using  a  coordinate-invariant  formulation,  where  problems  of  polyhedral 
inclusions  in  an  elastic  half-space  and  bi-materials  were  also  investigated.  In  addition,  specific  analytical  solutions  have 
been  obtained  for  polyhedral  inclusions  of  simple  shapes  such  as  cuboids  (e.g.,  Chiu,  1977;  Lee  and  Johnson,  1978;  Liu  and 
Wang,  2005)  and  pyramids  (e.g.,  Pearson  and  Faux,  2000;  Glas,  2001 ;  Nenashev  and  Dvurechenskil,  2010).  Also,  illustrative 
results  have  been  provided  for  dynamic  Eshelby  problems  of  cubic  and  triangularly  prismatic  inclusions  along  with 
spherical  and  ellipsoidal  ones  by  Wang  et  al.  (2005)  using  their  general  solution  for  the  dynamic  Eshelby  problem  for 
inclusions  of  various  shapes. 

However,  these  existing  studies  on  polyhedral  inclusion  problems  are  all  based  on  the  classical  elasticity  theory,  which  does 
not  contain  any  material  length  scale  parameter.  As  a  result,  the  Eshelby  tensors  obtained  in  these  studies  and  the  subsequent 
homogenization  methods  cannot  capture  the  inclusion  (particle)  size  effect  on  elastic  properties  exhibited  by  particle-matrix 
composites  (e.g.,  Vollenberg  and  Heikens,  1989;  Cho  et  al.,  2006;  Marcadon  et  al.,  2007).  Solutions  for  polyhedral  inclusion 
problems  are  also  important  for  describing  interpenetrating  phase  composites  reinforced  by  3-D  networks  (e.g.,  Poniznik  et  al., 
2008;  Jhaver  and  Tippur,  2009)  and  for  understanding  semiconductor  materials  buried  with  quantum  dots  that  are  typically 
polyhedral-shaped  (e.g.,  Kuvshinov,  2008;  Nenashev  and  Dvurechenskil,  2010).  These  materials  often  exhibit  microstructure- 
dependent  size  effects  whose  interpretation  requires  the  use  of  higher-order  continuum  theories. 

In  this  paper,  the  Eshelby-type  inclusion  problem  of  a  polyhedral  inclusion  prescribed  with  a  uniform  eigenstrain  and  a 
uniform  eigenstrain  gradient  and  embedded  in  an  infinite  homogeneous  isotropic  elastic  material  is  solved  using  a 
simplified  strain  gradient  elasticity  theory  (SSGET)  (e.g.,  Gao  and  Park,  2007),  which  contains  a  material  length  scale 
parameter  and  can  describe  size-dependent  elastic  deformations.  The  Eshelby  tensor  is  analytically  obtained  in  terms  of 
three  potential  functions,  two  of  which  are  the  same  as  the  ones  involved  in  the  counterpart  Eshelby  tensor  based  on 
classical  elasticity.  These  potential  functions,  as  three  volume  integrals  over  the  polyhedral  inclusion,  are  evaluated  by 
dividing  the  polyhedral  inclusion  domain  into  tetrahedral  duplexes.  Each  duplex  and  the  associated  local  coordinate 
system  are  constructed  using  a  procedure  similar  to  that  developed  by  Rodin  (1996)  based  on  the  algorithm  proposed  in 
Waldvogel  (1979).  Each  of  the  three  volume  integrals  is  first  transformed  to  a  surface  integral  by  applying  the  divergence 
theorem,  which  is  then  transformed  to  a  contour  (line)  integral  based  on  Stokes’  theorem  and  using  an  inverse  approach 
different  from  those  employed  in  the  existing  studies  for  evaluating  the  two  integrals  involved  in  the  classical  elasticity- 
based  Eshelby  tensor  for  a  polyhedral  inclusion. 

The  rest  of  this  paper  is  organized  as  follows.  In  Section  2,  the  general  form  of  the  Eshelby  tensor  for  a  3-D  arbitrary- 
shape  inclusion  based  on  the  SSGET  is  presented  in  terms  of  three  potential  functions  (volume  integrals).  The  expressions 
of  the  SSGET-based  Eshelby  tensor  for  a  polyhedral  inclusion  of  arbitrary  shape  are  analytically  derived  in  Section  3,  which 
is  separated  into  a  classical  part  and  a  gradient  part.  The  averaged  Eshelby  tensor  over  the  inclusion  volume  is  also 
analytically  evaluated  there.  Numerical  results  are  provided  in  Section  4  to  quantitatively  illustrate  the  position  and 
inclusion  size  dependence  of  the  newly  obtained  Eshelby  tensor  for  the  polyhedral  inclusion  problem.  The  paper  concludes 
in  Section  5. 

2.  Eshelby  tensor  based  on  the  SSGET 

The  SSGET  is  the  simplest  strain  gradient  elasticity  theory  evolving  from  Mindlin’s  pioneering  work.  It  is  also  known  as 
the  first  gradient  elasticity  theory  of  Helmholtz  type  and  the  dipolar  gradient  elasticity  theory  (e.g.,  Gao  and  Ma,  2010a). 
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According  to  the  SSGET  (Gao  and  Park,  2007;  Gao  and  Ma,  2010a),  the  Navier-like  displacement  equations  of  equilibrium  are 
given  by 

(X+ )i)uuj  +  fium-L2[().+ )i)uuj  +  )iulkk]  mm+fj  =  0  in  12,  (1) 

and  the  boundary  conditions  have  the  form: 

aHnj-^ijknk\j+^lijknkni\ini  =  c<  or  u>  =  u< 

HykOjOk  =  q,  or  uun,  =  g) 


on  8Q, 


(2) 


where  ).  and  /(  are  the  Lame  constants  in  classical  elasticity,  L  is  a  material  length  scale  parameter,  u,  are  the  components  of  the 
displacement  vector,])  are  the  components  of  the  body  force  vector  (force  per  unit  volume),  are  the  components  of  the  total 
stress,  <r  =cr, je.igiej,  f.iyk  are  the  components  of  the  double  stress,  n=/(,jke,-®ej(g>efo  and  t,-  and  q,  are,  respectively,  the  components 
of  the  Cauchy  traction  vector  and  double  stress  traction  vector.  Also,  in  Eqs.  (1)  and  (2),  12  is  the  region  occupied  by  the 
elastically  deformed  material,  512,  is  the  smooth  bounding  surface  of  12,  n,  is  the  outward  unit  normal  vector  on  512,  and  the 
overhead  bar  represents  the  prescribed  value.  In  addition, 


l^ijk  —  ^  ^ij.k  —  Hjik,  (Tjj  —  Tjj  fJ-ijkik  —  T,j  L  Tjjkk  —  Oji, 


(3) 


where  t,j  are  the  components  of  the  Cauchy  stress,  z  =  ryej<giej. 

When  the  strain  gradient  effect  is  not  considered  (i.e„  L= 0),  Hyk=  0  and  ery= t,j  (see  Eq.  (3)),  and  Eqs.  (1)  and  (2)  reduce 
to  the  governing  equations  and  the  boundary  conditions  in  terms  of  displacement  in  classical  elasticity  (e.g.,  Timoshenko 
and  Goodier,  1970;  Gao  and  Rowlands,  2000). 

Note  that  the  standard  index  notation,  together  with  the  Einstein  summation  convention,  is  used  in  Eqs.  (1 ),  (2)  and  (3) 
and  throughout  this  paper,  with  each  Latin  index  (subscript)  ranging  from  1  to  3  unless  otherwise  stated. 

Solving  Eq.  (1 ),  subject  to  the  boundary  conditions  of  u,  and  their  first,  second  and  third  derivatives  vanishing  at  infinity, 
gives  the  fundamental  solution  and  Green’s  function  based  on  the  SSGET.  By  using  3-D  Fourier  and  inverse  Fourier 
transforms,  the  fundamental  solution  of  Eq.  (1)  has  been  obtained  as  (Gao  and  Ma,  2009) 


rr  r  +oo 

uiW=  I]]  Gij(x-y)fj(y)dy, 


(4) 


where  x  is  the  position  vector  of  a  point  in  the  infinite  3-D  space,  y  is  the  integration  point,  and  Gy{  • )  is  Green’s  function 
(a  second-order  tensor)  given  by 

1  JA(x)dy-B(x)y],  (5) 


C«(x)  = 


167t/i(l-f) 


with 


A(x)  =  4(l-t>)i  (l-e  X/,l)>  B(x)  =  x+ 


2 Ll  2 Ll 


-x/L 


(6) 


When  L=0,  Eqs.  (5)  and  (6)  reduce  to  the  Green’s  function  for  3-D  problems  in  classical  elasticity  (e.g.,  Li  and  Wang,  2008).  In  Eqs. 
(5)  and  (6),  v  is  Poisson’s  ratio,  which  is  related  to  the  Lame  constants  X  and  through  (e.g.,  Timoshenko  and  Goodier,  1970) 
Ev  E 


k  = 


(1  +  v)(l—2v)  ’ 


f‘  = 


2(1  +  1')’ 


(7) 


where  E  is  Young’s  modulus. 

By  using  the  Green’s  function  method  entailing  Eqs.  (4),  (5)  and  (6),  the  general  expressions  of  the  Eshelby  tensor  based 
on  the  SSGET  can  then  be  obtained,  as  summarized  below. 

Consider  the  problem  of  a  3-D  inclusion  of  arbitrary  shape  embedded  in  an  infinite  homogenous  isotropic  elastic  body. 
The  inclusion  is  prescribed  with  a  uniform  eigenstrain  e*  and  a  uniform  eigenstrain  gradient  k*.  There  is  no  body  force  or 
any  other  external  force  acting  on  the  elastic  body.  The  disturbed  strain,  Sy,  induced  by  e*  and  k*  can  be  shown  to  be  (Gao 
and  Ma,  2009,  2010b) 

l-'ij  =  $ijkl‘'ki  +  Tijkhn  ,  (8) 

where  Syki  is  the  fourth-order  Eshelby  tensor  having  36  independent  components,  and  Tykim  is  a  fifth-order  Eshelby-like 
tensor  with  108  independent  components.  Eq.  (8)  shows  that  e  (  =  £ye,igiej)  is  solely  linked  to  e*  in  the  absence  of  k*  (i.e„ 
the  classical  case)  and  is  fully  related  to  k*  when  e*=0.  The  fourth-order  Eshelby  tensor  has  been  obtained  as 


Sijkl  —  Sijkt+$jkl’ 


_ 
kl  — 


1 


871(1  —V) 


cC  _  1 

87t(i  -u) 


[<&,ijkl—2vAijdki—('l -v)(Ajj5ik+  A  kjdu  +  A  ,idjk  +  A  kiSji)], 


[2  ofySiti+o -v)(f  jid  ik + r  udjk + r  jk8u + r  ik5ji) + 2i?(Atijki- r  pO], 


(9a) 


(9b) 


(9c) 
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where  Sfjkl  is  the  classical  part,  Sykl  is  the  gradient  part,  <5y  is  the  Kronecker  delta,  and 


<Z>(x)=  /  jx— y|dy, 

(10a) 

JQ 

MX)-  £iv  vidy- 

(10b) 

e-|x-y|/I 


Fix)  -  dy 

(10c) 

JQ  x-y| 

are  three  scalar-valued  potential  functions  that  can  be  obtained  analytically  or  numerically  by  evaluating  the  volume 
integrals  over  the  domain  Q  occupied  by  the  inclusion,  with  jxj  =x=(xkxkY12  and  y  (e(2)  being  the  integration  variable. 
Note  that  the  first  two  potential  functions  given  in  Eqs.  (10a)  and  (10b)  are  the  same  as  the  ones  involved  in  the  Eshelby 
tensor  based  on  classical  elasticity  (e.g.,  Mura,  1987;  Nemat-Nasser  and  Hori,  1999;  Li  and  Wang,  2008),  whereas  the  third 
one  defined  in  Eq.  (10c)  results  from  the  use  of  the  SSGET.  It  should  be  mentioned  that  <P(x)  in  Eq.  (10a)  and  yl(x)  in 
Eq.  (10b)  are,  respectively,  known  to  be  a  biharmonic  potential  and  a  Newtonian  potential  (e.g.,  Li  and  Wang,  2008),  while 
a  variant  of  F(x)  in  Eq.  (10c)  is  called  the  Yukawa  potential  in  physics  (e.g.,  Rowlinson,  1989). 

Eqs.  (10a),  (10b)  and  (10c)  show  that  among  the  three  potential  functions,  only  the  third  one,  F(x),  involves  the  length 
scale  parameter  L.  It  then  follows  from  Eqs.  (9c),  (10b)  and  (10c)  that  Sykl  depends  on  L,  while  Sykl,  expressed  in  terms  of 
A(x)  and  <£(x)  only  according  to  Eqs.  (9b),  (10a)  and  (10b),  is  independent  of  L.  Also,  it  is  seen  from  Eqs.  (9c),  (10b)  and 
(10c)  that  Sfjkl  =  0  when  L= 0  (and  thus  F(x)= 0),  thereby  giving  Sp;  =  S^kl  (from  Eq.  (9a)).  That  is,  the  Eshelby  tensor  derived 
using  the  SSGET  reduces  to  that  based  on  classical  elasticity  when  the  strain  gradient  effect  is  not  considered. 

The  fifth-order  Esheiby-like  tensor  T,  which  relates  the  eigenstrain  gradient  k*  to  the  disturbed  strain  £  in  the  elastic 
body  (see  Eq.  (8)),  can  be  shown  to  be 

]} 

Tijklm  =  -j  _ ^  [2vYf,ijm(5w  +  (l  ~ P)i  ,Fj,llj(')jk  -f  Jrnj(\k  +  ll\kj,ui)ji  -f  *F ikmj^il)~FI (H) 

where 

<P(x)  =  A-F,  n(x)  =  <F+2L2(A-D,  (12a,  b) 

with  the  scalar-valued  potential  functions  t>{x),  A(x)  and  f(x)  defined  in  Eqs.  (10a),  (10b)  and  (10c).  It  is  clear  from 
Eq.  (11 )  that  T  vanishes  when  L=  0.  Then,  with  Sykl  =  Sykl  as  discussed  above,  Eq.  (8)  simply  becomes  e,j-  =  Sykle^  when  L= 0, 
which  is  the  defining  relation  for  the  Eshelby  tensor  based  on  classical  elasticity  (Eshelby,  1957),  as  expected. 

Eqs.  (9a),  (9b),  (9c)  and  (11)  provide  the  general  formulas  for  determining  S yW  (  =  Sfjkl+Sfjkl)  and  Tykim  for  an  inclusion  of 
arbitrary  shape  in  terms  of  the  potential  functions  /l(x),  <P(x)  and  F(x)  defined  in  Eqs.  (10a),  (10b)  and  (10c).  For  the  cases 
of  a  spherical  inclusion,  a  cylindrical  inclusion  and  an  ellipsoidal  inclusion  in  an  infinite  elastic  medium,  analytical 
expressions  have  been  obtained  for  /l(x),  <F(x)  and  T(x)  and  thus  for  the  Eshelby  tensor  (Gao  and  Ma,  2009,  2010b;  Ma  and 
Gao,  2010).  The  more  complex  case  of  a  polyhedral  inclusion  of  arbitrary  shape,  for  which  /l(x),  <P(x)  and  f(x)  are  difficult 
to  evaluate  analytically,  is  examined  in  this  study. 

3.  Polyhedral  inclusion 

The  problem  of  an  arbitrary-shape  polyhedral  inclusion  in  an  infinite  elastic  body  has  been  analytically  studied  by  Rodin 
(1996),  Nozaki  and  Taya  (2001)  and  Kuvshinov  (2008)  using  classical  elasticity.  The  Eshelby  tensor  for  this  problem  is  derived 
here  using  the  SSGET-based  general  formulas  and  a  new  method  for  evaluating  the  potential  functions  /l(x),  <P(x)  and  T(x). 

3.1.  Eshelby  tensor 

Consider  an  arbitrarily  shaped  polyhedral  inclusion  embedded  in  an  infinite  homogeneous  isotropic  elastic  material. 
The  polyhedral  inclusion  has  p  faces  and  is  prescribed  with  a  uniform  eigenstrain  £*  and  a  uniform  eigenstrain  gradient  k*. 

The  p-faced  polyhedral  domain  occupied  by  the  inclusion  can  be  divided  into  tetrahedral  duplexes  originated  from  a  chosen 
(arbitrary)  point  x  (Waldvogel,  1979;  Rodin,  1996).  Each  duplex  can  be  further  divided  into  two  simplexes,  each  of  which  is  a 
tetrahedron  with  three  of  its  four  faces  being  right  triangles  (see  Fig.  1 ).  The  four  vertices  of  each  of  the  duplexes  are,  respectively, 
the  projection  point  of  x  on  a  polyhedral  surface  (i.e.,  x(),  two  adjacent  vertices  on  this  surface  (i.e.,  V(|  and  Vjj ),  and  the  point  x 
itself.  For  each  of  these  duplexes,  a  local  Cartesian  coordinate  system  is  constructed,  with  point  x  being  set  as  the  origin.  The  three 
orthogonal  axes  of  the  local  coordinate  system  are  denoted  by  X,  p  and  (,  respectively.  The  coordinates  of  the  two  vertices 
Vjl  and  Vjj  on  the  Jth  edge  of  the  /th  surface  are,  respectively,  given  by  (bp,  Ip ,  a;)  and  (bp,lp,a{),  as  shown  in  Fig.  2. 

In  Fig.  2,  Xp,  njj  and  ^;°  are  the  unit  vectors  associated  with  the  local  coordinates  Xji,  i)ji  and  flt  y  is  an  arbitrary  point  on 
the  Jth  edge  of  the  /th  surface,  r  is  the  position  vector  of  y  relative  to  the  origin  x  (i.e.,  r=y-x),  and  ij  is  the  projection  of 
r  on  the  /th  surface.  The  usual  Cartesian  coordinates  (x-i,  x2,  x3)  are  used  in  the  global  coordinate  system  having  (elt  e2,  e3) 
as  the  associated  base  vectors. 
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a  b 


Fig.  1.  A  polyhedron  represented  by  duplexes:  (a)  a  polyhedron  (with  five  duplexes  shown);  (b)  a  duplex  and  the  associated  local  coordinate  system 
constructed  from  an  arbitrary  point  x. 


To  obtain  the  Esheiby  tensor  for  the  polyhedral  inclusion  using  Eqs.  (9a),  (9b)  and  (9c),  the  three  potential  functions 
<£(x),  yl(x)  and  F(x)  defined  in  Eqs.  (10a),  (10b)  and  (10c)  are  first  evaluated  over  the  polyhedral  domain  Q  using  an 
approach  different  from  those  employed  in  Rodin  (1996),  Nozaki  and  Taya  (2001)  and  Kuvshinov  (2008)  for  evaluating 
<P(x)  and  zt(x)  involved  in  the  classical  elasticity-based  Esheiby  tensor,  as  shown  next. 

For  a  sufficiently  smooth  function  M(x-y),  the  use  of  the  divergence  theorem  gives 

where  p  is  the  number  of  surfaces  of  the  polyhedron,  and  (C,)fc  is  the  kth  component  of  the  unit  outward  normal  vector  on 
the  /th  surface  8Qh 

To  transform  the  surface  integral  in  Eq.  (13)  to  a  contour  (line)  integral,  let 

M  =  (Vxm)'(f,  (14) 


where  m  is  a  yet-unknown  vector  located  on  the  /th  surface  of  the  polyhedron,  and  V  x  m  denotes  the  curl  of  m.  Using  the 
Stokes  theorem  then  yields,  upon  applying  Eq.  (14), 


//, 


MdS(y)  = 

dC2/  J  =  1  ' 


m  ■  rjjjdl, 


(15) 


where  rjjJ  is  the  unit  vector  along  the  Jth  boundary  edge  Cj,  of  the  /th  surface. 
Now,  write 


m  =  tfx 


if 

-kg(r) 


(16) 
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where  rf(  =  ^Jr2—aj  =  | rf  | )  is  the  length  of  the  projection  of  r  on  the  /th  surface,  and  g(r)  is  a  function  of  r  (=  jr|)  yet  to  be 
determined.  Substituting  Eq.  (16)  into  Eq.  (15)  leads  to 

[f  MdS(y)  =J'b]l  [  ^-dl,  (17) 

JJ  sa,  fr\  JCji  r/ 

where  bjt  =  rf  ■  is  the  distance  from  point  X;  (the  projection  of  point  x  on  the  /th  surface)  to  the Jth  edge  Cj /  (see  Fig.  2).  For 
each  specific  function  M,  a  different  expression  of  g(r )  and  thus  m  can  be  determined,  as  shown  next  for  the  three  cases 
representing  the  integrands  of  the  potential  functions  $(x),  /l(x)  and  T(x)  defined  in  Eqs.  (10a),  (10b)  and  (10c). 

For  M=r=  |y— x|  (corresponding  to  <P(x)),  Eqs.  (14)  and  (16)  gives 


r  =  g(r)[  V- J 


+  [Vg(r)]-  ' , 

ri 


(18) 


where  V  is  the  gradient  operator,  and  use  has  been  made  of  the  identity:  a  x  b  x  c=(a  ■  c)b  —  (a  ■  b)c,  with  a,  b,  c  being 
arbitrary  vectors  and  “  x  ”,  “  •  ”  representing  the  cross,  dot  products,  respectively.  After  carrying  out  the  differentiation  and 
dot  product  operations,  Eq.  (18)  can  be  further  simplified  to 

r[ 

r 


r-f-Httr ,4. 


(19) 


where  g,(  =  dg/dr)  is  the  first  derivative  of  g  with  respect  to  r.  The  solution  of  Eq.  (19)  reads 
r3-a? 


g4>(r)  = 


^-n3 
3i?  ’ 


(20) 


where  a7  (  =  r ■£, )  is  the  distance  from  point  x  to  the  /th  surface  (see  Fig.  2),  and  g&{r)  denotes  the  function  g(r)  for  the  case 
with  M=r. 

Similarly,  it  can  be  shown  that 


ga(r)  = 


_A_ 

r+a, 


when  M=l/r=l/|y-x|  (corresponding  to  /l(x)),  and 
L(e-°  i/L~e-r/L) 


grin  = 


(21) 


(22) 


when  JV/  =  e_r''1/r  =  e“ly_x[/i/|y-xj  (corresponding  to  T(x)). 

Using  Eqs.  (13),  (17),  (20),  (21)  and  (22)  in  Eqs.  (10a),  (10b)  and  (10c)  then  leads  to,  with  the  local  coordinate  axis  t] 
being  along  the  Jth  edge, 


V  4  rl, 

A,= -££(£, W 

J=1J  =  1  Jlj> 


'/  (af  +  bJ,+if)3/2-af 


di /, 


A,  =  - 


3  (bf,  +  t/2) 

EE«?*/r  ;  <*!■ 

f  =  i/  =  i  Jv  Jaf  +  bj,  +rj2+a, 

(J  ]_(e~a,/L— e_\/a'  +bi‘ +l*2/1) 


bj,  +  ri2 


dij, 


(23) 


(24) 


(25) 


where  lj[  and  lj,  are,  respectively,  the  coordinates  of  the  two  vertices  V/T  and  v/7  on  the  Jth  edge,  with  !,}  being  positive  and 
IJf  negative  (see  Fig.  2). 

The  integrals  in  Eqs.  (23)  and  (24)  can  be  exactly  evaluated  by  direct  integration  to  obtain  the  following  closed-form 
expressions: 


*4  =  -  IE  (tf).  \  ^  s/af+b2+(ljl)2+  i 


tan 


a, Ip 


a. 


—  h-tan  hr-  + 


i  AT  \  3  ajbji  +  bp 


In 


b„y/af  +  bj,  +  (lj T)2 

'jT  +  ^a}  +  bji+0ji)2 


\Ja2+b 
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_  "it  ^a'2  +  +<'l]i)2~  T 


tan 


3 

p  i 


Of  J// 


fa;Va'2 +bji+('lji^2 


+!-©-3^ln 


lj,  +  ^/a2  +  bjI  +  (lj,)2 


aj  +  bj, 


—  EE  =-EE  (tfw(<)+- -(<)“]• 

/=!_/=!  /=!_/=! 


Ai-  =  -  E  a,tan  1 

i=ij=i  i 


a/'jt 


bJi  \J  °2  +  bp  +  Op  >2 


+bj,  In 


+  \/a?  +  bji+^n  *2 


-a(tan 


-«/tan  1  f  ] 
bJi) 


b„Jaj  +  bj,+(ljif 


+  a;tan  1  (  )  — £>//  In 

W 


lji  +  \/a}  +  tfi  +  Oji)2 


-EE(M  =  -EE(*(  4»+-(4n. 

/  =  l  j  =  l  /  =  l  j  =  l 

where  (</>'') +  ,  (d^'r,  As{,  (/tf  j+and  (/tJ/r  are  functions  defined  by 

4  =  d^iat.bjhlp  Op)  =  (d^ )+(aI,bJI,Ijl )-<<P,Jr<a,,bJI,lp), 


<«?>  *-‘»j?V“f  +  6?.+('J)2  +  T 

3afbjI+bj, 


tan 


ai(H 


+  - 


In 


3  1  bji \jal +bp+0p)2 

Ijj  +  l/0?  +t>;/  +  ((//")' 


-ftan-Wf 
3  l  bj, . 


\[tf  +  b] 


^ +  bf' + (£)2  +  | 

'icifbji+bj. 


tan 


bj,Jaj  +  bj,+(.lj,)2 


+  - 


3  In  ^  +  +  ^  +  ^Jl  *2 


6  ya/2  +  faj/ 

4  =  4WMi -1/7)  =  (4)+(a/.fa;/./;T  )-(4)-(a,,h;,,i;7). 


(AJ-{) +  =  a(  tan  1 
(4)-  =  Q/tan-1 


aQp 

bJiy/tf  +  bp+0p)2 


_i  //,t\ 

-a(tan  *2-  +  fay/  In 
\bJiJ 


Ip  +  l/a;2  +bf,  +  (lp  )2 


aj  +  bj( 

(//  +  \Z<32  +  b2/  +  (/;/  )2 


app 

[tyf^/a?  +b2/  +  (/;;) 

Note  that  e~r/i  can  be  written  as  a  power  series 

£  W- 

n  =  0 

Using  Eq.  (29)  in  Eq.  (25)  then  leads  to 


,/  I  ,,  , 
a(tan  Ur-  +  b(,ln 

V'V 


\[tf+b] 


o=-EI(At|^ 

j  =  u  =  i 


tan  1 1  -tan  1  (E  I 


v^(-i)n+,!?[flf+t$+(!p)2]n/2 
+  „  =  o  ^  ^ 

_  ^  (-l)"+1/j7[fl?+bj,+(ij7)2]n/2 

b"n!  bp 


fi+)2  ]“n/2 

1  +  ^ — ,2 

ai  +  bji. 


n  =  0 
P  Q 


1  + 


Op)2  ' 

tf+tfl. 


-n/2 


1  n  3  (Ip)2  Op)2 

r  r  'v  a2+b2 •  h2 

1  _n  3  (fj7)2  (|;7)2 

[2’  2’  '2'  a2  +  b2’  b2  J 


—  EE  = -  E  E  (c?)i[(4)+-(4)"]. 

/=!/=!  /  =  1 J  =  1 


(26) 


(27) 


(28) 


(29) 


(30) 
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where  r1',  (rJ1,)+  and  (F^)  are  functions  defined  by 

rJ!  =  ,/j7)  =  (/’{')+ (^  V'jT  )-(rJ/r(an  V'/7). 

(— l)n+1  lji[aj  +bji+(lji)2]n/2 


(r^=L^  tan-^)+E^ 


1  n  3  (IjlV 

r~r  'T  aj+by 


i  + 


(/;T)2rn/2 


af  +  tfL 


djjr 


UJI 

(-1)"+1  ljI[af  +  tfl+(lJI)2]n/2 


UH 


1  + 


('„> 


2  n-n/2 


a?  +  tf, 


F, 


1  _n  1  3  _  (Z,/)2 


(/,7)2 


af  +  bj ,  bj, 


and  Fi  is  the  first  Appell  hypergeometric  function  of  two  variables  given  by 

Fi(a,b,c,d;x,y)  =  £  (°)"+n(^m,(n?nxn,yn. 

m,n  =  0  Wm  +  n"1-’1- 


(31) 


(32) 


with  C/)m  being  the  Pochhammer  symbol  representing  the  following  rising  factorial: 

(f)m=ri^^=f(f+\)(f+2)---(f+m-\).  (33) 

Note  that  Eqs.  (26),  (27)  and  (30)  are  applicable  to  both  the  interior  case  with  x  being  inside  the  polyhedral  inclusion 
(i.e.,  xefl)  and  the  exterior  case  with  x  being  outside  the  inclusion  (i.e.,  x$Q).  For  the  former  a,  is  a  positive  value,  while  for 
the  latter  a;  is  a  negative  value.  The  similarity  and  difference  identified  here  between  the  interior  and  exterior  cases  can  be 
seen  from  Fig.  3,  where  how  the  duplex  in  each  case  is  constructed  is  schematically  shown. 

It  should  be  mentioned  that  no  attempt  is  made  here  to  obtain  the  expressions  of  the  potential  functions  <P(x),  /L(x)  and 
F(x)  from  Eqs.  (26),  (27)  and  (30),  since  only  the  second  and/or  fourth  derivatives  of  these  functions  are  involved  in  the 
general  expressions  of  the  Eshelby  tensor  given  in  Eqs.  (9a),  (9b)  and  (9c). 

Note  that  for  a  smooth  function  F(x)  =  F(a(,  bjiJjlJji)  =  F+(ai,bjiJji)—F  (a/.bji,  5)  the  use  of  chain  rule  gives 

8 F  8F  da i  dF  dbji  dF+  dlj [  8F~  8ljt 
■'  ~  dXj  ~  da,  dXj  +  dbj,  dxt  +  dip  8xt  8ljt  dXj  ’ 

where  the  parameters  a,,  bp.  Ip  ,  IJj  are  related  to  x  through 

a,  =  (v^-xk)(C°)k, 


(34) 

(35a) 


bji  =  (vk  -xk)(Xp)k, 


(35b) 


l/T  —  (vk  ~xk)(t]p)k. 


(35c) 


Fig.  3.  Duplex  and  parameters  a/t  bji,  Ijl,  lj,  for  (a)  xeQ  and  (b)  x$Q. 
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hi  =  (pk  ~ *k)(riji)k . 


(35d) 


where  xk,  and  v k  are,  respectively,  the  coordinates  of  the  points  x,  V,}  and  v;7  in  the  global  coordinate  system,  and  (£/)t, 

(/,p)k  and  are  the  components  of  the  unit  base  vectors  7.°,  and  t|jj  in  the  global  coordinate  system. 

It  then  follows  from  Eqs.  (34),  (35a),  (35b),  (35c)  and  (35d)  that 


_  SF  o  8F  o,  ( SF+  8F  >  0 


(36a) 


&F,r  o 


d3F, 


F  m  =  - 


d3F+  c?F- 


c?F 


■0\  /  t  0  /  ^0-. 


a(/,t  r  5(/;7)JJ 


/  i  o 


~  d3F+  c?F- 


+  (£;  )i(Cj  )j(^7j/)fe] — : 


da,db ’], 

+  (^ji)i(rlji)j(^ji)k+(Aji)i(Aji)j(t1p)k\— 

a3F+  a3F 


da}  dip  da}  dip  J 


K^¥C?)A+(C?)KA(tf)t 

53F+  ^F' 


ebpipy 


53F+ 


^F- 


[(£/  )iOljl)j(rfl)k  "F  (*7//)i(£/  Ijibjlh  +  i'ljllMjlljiCl  )kl 


o-i  /r°'ivviO\ 


1,0  I  ,  /rO, 


aa,a(/+)2  aa,a(/;7)2J 

[(Z/j)i(t?/j  )j{>lp)k  +  (>ljl  hi^jllji'ljl  Ik  +  0ljl  liOljl  )j  (Aji.  )fc] 


dbpd(lj})2  dbpd(lp)2 

_  (, dafdbpdlj 7  “  dafdbpdl j)  K^WC/%07,% 

+  (^// )i(t/;/)j(C/  )k  +  (t/jj  ),(£/  )j(^ji)k  +  (>lji)i(^]i)j(Ci  )fc]- 


(36b) 


Using  Eqs.  (26),  (27),  (30),  (36a)  and  (36b)  in  Eqs.  (9b)  and  (9c)  will  lead  to  the  final  expressions  of  the  Eshelby  tensor 
for  the  p-faced  polyhedral  inclusion  as 

i  p  ? 

cC  -  1  ' 

®w“  871(1 


+  (Sf)J/[(C?)fc(C?)J<5„ +(C?)i(C?)j<5lfc] + (£)/£?),($/«  +  (Sf)Jf  [(C?)f(A°)fc<5j7 +(C?),a°)l«5Jfc] 

+  (^6  )y/[(7(//)fc(C/  )j^a + (^ji)i(C/  )j<5ifc] 4- (5y )j/(C/  ([(t/J/lj^M  +  fSglj/KCj  )i(^7//)fe<5ji + (C/  IMjilibjkl 

+  (Sg);/[0/j})|;(£j  )j<5il +(>?;/)/(£/ )j<5»]+(S^o)j/[(£j  )k(C?)l  +  (£/  )i(£?)j(^/()fc(£f  )l+(£f  )i(C/  )j(£/  M^y/)/] 

+  (^  )y/[(Cy0)i(»7J0/)J(C/°)k(C?),  +  )k(C?),  +  (C/0)^^0)/^0^)^0/  >/]  +  (^2 


£,V°)y(^ 


cC  VKrf)i(tf  ),t 


r°)t(</)J(C?J 
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+  (£/  )i(Zy/)y(Cj  )k07y/)|  +  (Cj  )i(tIjl)j(£l)k(^jih+(C<l)iOljl)j(tfl)k(Cl)l])’  (37a) 


where 
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(37b) 
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tfkl  =  8na -V)i  + 


+  (^4  )ji(Ci  )i(^]i)jSkl+(Sl^)jl[{Cl  )i(lji)kSjl  +  (fj  +  (S6)j/[(^)([(C/  l/^i! +  ('$)((£/ 

+  (^7 )j/(C/  )iOlji)jdkl  +  (Sg  );/[(C;  )j(bji  )k^ji I  +  (Cl  )i(,1jl)l^jk]  +  (Sg  )/;[(>/;;  )fc(C/  )j^il  +  Olji )f(Ci  )j^i(r] 

+ (Sfo  + (CfiKC,0)^^  ^(G0), + (C?)i(C/°)J(C5,)fc(  >,] + (Sf,  ^[(C/0)^^0/  )J(C/°)fc(C/0)/ 
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+  (Cl  )i(rl] | +  (Cl  ¥'7/I )j(J7jl)A:(-C/i)l]  +  (S^5  )ji[(Ci  )i(>1ji )j(^JI hi^JI  h  +  (Cl  )i('^/l)j(^7//  )lc(^//)l  +  (Cl  )i(-C/i)j(^/i  )k(b;i ) 
+  (^?6);i[(Cl  )i(Ci  )j(^-JI^fc(f7j/)/  +  (Ci  )i(Ci  )j(f7jl)lc('C/i)( +  (Ci  )i(-lj})j(>7j()fc(Cj  )l  +  (Ci  )i(2//)j(Ci  )k(Tl°l) 
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In  Eqs.  (37b)  and  (37d),  (<Ply+ ,  (<b,/)“,  /l^1,  (/t-,/)  +  ,  F]j,  (r]l)+  and  (F-^r  are  defined  in  Eqs.  (28)  and  (31). 

It  should  be  mentioned  that  the  classical  part  Sp,  in  Eqs.  (37a)  and  (37b)  depends  only  on  Poisson’s  ratio  v  and  cannot 
account  for  the  inclusion  size  effect,  noting  that  cb^,  (<b,/)+,  (<b^r,  /if,  (/f^)*  and  (/1J|,)~  involved  in  Eq.  (37b)  do  not  contain 
the  material  length  scale  parameter  L  (see  Eq.  (28)).  However,  the  gradient  part  Sykl  in  Eqs.  (37c)  and  (37d)  can  capture  the 
inclusion  size  effect,  since  Eqs.  (37c)  and  (37d)  as  well  as  the  expressions  of  rJ(,  (rJ1,)  +  and  (F^r  (see  Eq.  (31))  contain  the 
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parameter  L  in  addition  to  Poisson’s  ratio  v.  Clearly,  when  L= 0  (i.e.,  in  the  absence  of  the  strain  gradient  effect),  Sjjk/  =  0 
according  to  Eqs.  (37c),  (37d)  and  (31 ),  thereby  resulting  in  Sijkl  =  S^kl  from  Eq.  (9a).  That  is,  the  SSGET-based  Eshelby  tensor 
reduces  to  its  counterpart  based  on  classical  elasticity  when  the  strain  gradient  effect  is  not  considered.  Also,  it  is  seen  that 
the  expressions  of  the  classical  elasticity-based  Eshelby  tensor  in  Eqs.  (37a),  (37b)  and  (28)  derived  here  are  more  compact 
than  those  given  in  Nozaki  and  Taya  (2001). 

The  expressions  of  the  Eshelby  tensor  Syu  in  Eqs.  (9a),  (37a),  (37b),  (37c)  and  (37d)  are  derived  for  a  p-faced  polyhedral 
inclusion  of  arbitrary  shape.  For  simple-shape  inclusions,  more  explicit  expressions  can  be  obtained  for  Syki. 


3.2.  Averaged  Eshelby  tensor 


The  volume  average  of  the  position-dependent  Eshelby  tensor,  Sp(,  is  given  by 
Sijkl  =  X!  1212  f[f  ,lj,)dV, 

v^M=1N=\I=1J=r  JJQnm 


(38) 


where  (Sjvm)ijW  is  the  Eshelby  tensor  at  point  x  inside  QNM  presented  in  Eqs.  (9a),  (37a),  (37b),  (37c)  and  (37d),  Va  is  the 
volume  of  the  polyhedral  inclusion  Q,  QNM  is  the  region  occupied  by  the  duplex  formed  by  the  origin  (point  O)  of  the  global 
coordinate  system,  the  projection  of  point  O  onto  the  Mth  polygonal  surface  (i.e.,  Om)  and  two  vertices  on  the  Nth  edge  of 
the  Mth  surface  (i.e.,  \I^M  and  V^M),  and  n  is  the  number  of  edges  on  the  Mth  surface.  Note  that  this  duplex  QNM  is  different 
from  that  formed  by  point  x,  its  projection  onto  the  1th  polygonal  surface  (i.e.,  x,)  and  two  vertices  on  the  Jth  edge  of  the  1th 
surface  (i.e.,  V,}  and  Vjj),  as  shown  in  Fig.  4. 

For  the  NMth  duplex  QNM  originated  from  point  O,  the  local  Cartesian  coordinate  system  ( xNM ,  >Inm,  Cm)  can  be  chosen  in  a 
way  similar  to  what  was  done  earlier  (see  Fig.  2).  Then,  the  coordinates  of  the  vertices  of  the  duplex  on  the  Jth  edge  of  the  Jth 
surface  and  of  an  arbitrary  point  x  within  the  NMth  duplex  QNM  in  the  (2NM,  i)NM,  Cm)  local  coordinate  system  can  be  identified 
as  (v\INM+ ,  pfM+,  vfM+),(<NM-,  vfM-)and(xTM,  xf™,  xfM),  respectively.  Also,  the  base  vectors  n°  and 

of  the  local  coordinate  system  attached  to  the  duplex  Q}1  originated  at  x  can  be  expressed  in  terms  of  the  base  vectors  r$M 
and  Cm-  It  then  follows  that  the  parameters  for  the  duplex  Qj,  can  be  determined  as 

a,  =  (<NM+  (39a) 

bJt  =  (v>lNM+  ~^NM)(2°)NkM  =  (i/k'NM--x''NM)a°)r.  (39b) 


Fig.  4.  Duplexes  and  the  corresponding  local  coordinate  systems  constructed  from  an  arbitrary  point  x  and  from  the  origin  O  of  the  global  coordinate 
system,  respectively. 
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£ =(v,r+-<NMHfkM. 

(39c) 

ij,=(v>r--xTH)NkM , 

(39d) 

where  (Aj,)kM,  (i] and  represent,  respectively,  the  fcth  components  of  the  unit  vectors  Xp,  Tjjj  and  £,°  in  the  local 

coordinate  system  {XNM,  >1nm,  Cm)  with  the  base  vectors  X°NM,  and  Cm- 

Using  Eqs.  (39a),  (39b),  (39c)  and  (39d)  in  Eq.  (38)  yields 

_  ^  P  ^  n  ^  P  ^  Q  ^  r  r  r 

S, jkl  =  w—  ^2  NM-Ciw).  1)1  (/'-NM,'/nM'Cm).  Ijl  UnM.'I  ■  (40) 

V 12  M  =  1  N  =  1  1  =  1 J  = 

This  general  formula  can  be  used  for  a  polyhedral  inclusion  of  arbitrary  shape. 

For  a  polyhedral  inclusion  that  is  symmetric  about  the  global  coordinate  axes  x3,  x2  and  x3,  only  one  eighth  of  the 
inclusion  needs  to  be  considered  and  the  global  coordinate  system  can  be  used  in  all  computations.  The  one-eighth 
polyhedral  domain  can  be  divided  into  several  sub-polyhedra  with  their  top  and  bottom  surfaces  parallel  to  the  XiX2-plane, 
and  the  volume  integral  over  each  sub-domain  can  be  evaluated  by  direct  integration  using  the  global  coordinate  system. 
Also,  only  the  global  coordinates  of  all  vertices  need  to  be  determined,  and  the  unit  vectors  qjj  and  C;°  in  the  local 
coordinate  system  can  be  expressed  in  terms  of  the  base  vectors  (i.e.,  elt  e2,  e3)  in  the  global  coordinate  system.  As  a  result, 
Eq.  (40)  can  be  simplified  to 

8  ^  P  ^  Q  f  f  f 

Sijkt  (ST),jH[a,(x1,X2,x3),  bj)(x,,x  2,x3),  /,}  (XuX2,x3),  /)7(xi  ,x2  ,x3  )]dX)  dx2  dx3 ,  (41) 

*  ®  T  =  1 1  =  1J  =  ]OjJqt 

where  t  is  the  number  of  sub-polyhedra  in  the  one  eighth  of  the  polyhedral  inclusion,  and  ( Sj)ijki  is  the  Eshelby  tensor  at 
point  x  inside  Qr  given  in  Eqs.  (9a),  (37a),  (37b),  (37c)  and  (37d). 

4.  Numerical  results 

To  illustrate  the  general  formulas  of  the  Eshelby  tensor  for  a  p-faced  polyhedral  inclusion  of  arbitrary  shape  derived  in 
Section  3,  three  types  of  polyhedral  inclusions  (i.e.,  cubic,  octahedral  and  tetrakaidecahedral)  shown  in  Fig.  5  are 
quantitatively  studied  in  this  section.  Cuboids  are  the  first  polyhedral  inclusions  investigated  using  classical  elasticity  (e.g., 
Chiu,  1977;  Lee  and  Johnson,  1978;  Liu  and  Wang,  2005).  A  tetrakaidecahedron  can  be  generated  by  uniformly  truncating 
the  six  corners  of  an  octahedron  and  is  known  to  be  the  only  polyhedron  that  can  pack  with  identical  units  to  fill  space  and 
nearly  minimize  the  surface  energy  (e.g.,  Li  et  al.,  2003).  Tetrakaidecahedral  cells  have  been  frequently  used  to  represent 
foamed  materials  and  interpenetrating  phase  composites  (e.g.,  Li  et  al„  2003,  2006;  Jhaver  and  Tippur,  2009). 

Two  components,  S1311  and  S1212,  of  the  Eshelby  tensor  at  any  point  x  inside  each  polyhedral  inclusion  with  various 
sizes  are  evaluated  using  Eqs.  (9a),  (37a),  (37b),  (37c),  (37d),  (28)  and  (31)  and  plotted  to  demonstrate  how  the 
components  change  with  the  position  and  inclusion  size.  Also,  how  the  average  Eshelby  tensor  component  Sun  varies 
with  the  inclusion  size  is  presented  here,  which  is  computed  using  Eq.  (41).  For  illustration  purposes,  Poisson’s  ratio  v  is 
taken  to  be  0.3  and  the  material  length  scale  parameter  L  to  be  17.6  pm  in  the  current  numerical  analysis,  as  was  done 
earlier  (Gao  and  Ma,  2009,  2010a,b;  Ma  and  Gao,  2010,  2011). 

The  distributions  of  Smi  for  the  cubic,  octahedral,  and  tetrakaidecahedral  inclusions  along  the  x3  axis  predicted  by  the 
current  model  are  shown  in  Figs.  6-8,  where  the  values  of  S^in  are  also  displayed  for  comparison. 

It  can  be  seen  from  Figs.  6-8  that  the  classical  part  S,ni  (based  on  classical  elasticity)  varies  with  the  position  of  x 
within  each  polyhedral  inclusion  rather  than  uniform,  which  shows  that  the  Eshelby  conjecture  is  true  for  the  three 
polyhedral  inclusion  shapes  considered  here.  Also,  it  is  found  that  for  each  of  the  three  inclusion  shapes  S3111  at  a  given 
value  of  X)/R  is  the  same  for  all  values  of  R/L,  confirming  the  inclusion  size-independence  of  the  classical  part  of  the 
Eshelby  tensor,  which  is  noted  near  the  end  of  Section  3.1.  In  addition,  for  all  three  polyhedral  inclusion  shapes  considered, 
it  is  observed  from  Figs.  6-8  that  when  the  characteristic  inclusion  size  R  (see  Fig.  5)  is  small  (compared  to  the  length  scale 


Fig.  5.  Three  types  of  polyhedral  inclusions:  (a)  cubic,  (b)  octahedral  and  (c)  tetrakaidecahedral. 
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Fig.  6.  Variation  of  S1H1  along  thex!  axis  inside  the  cubic  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  edge  length  (see  Fig.  5(a)). 
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Fig.  7.  Variation  of  Sim  along  the  Xi  axis  inside  the  octahedral  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  edge  length 
(see  Fig.  5(b)). 
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Fig.  8.  Variation  of  Sim  along  the  Xi  axis  inside  the  tetrakaidecahedral  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  cell  height 
(see  Fig.  5(c)). 
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Fig.  9.  Variation  of  S1212  along  the  x-i  axis  inside  the  cubic  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  edge  length  (see  Fig.  5(a)). 


parameter  L,  e.g.,  R/L= 2),  the  strain  gradient  part  which  is  the  difference  between  Sun  and  Sj111  (i.e.,  S, , , ,  = 
Sn  n  -Sn  j  j )  and  is  displayed  as  the  vertical  distance  between  the  S'),  t  j  curve  and  each  St  ]  11  curve  in  Figs.  6-8,  is  significant 
and  should  not  be  neglected.  However,  as  the  inclusion  size  becomes  larger,  the  values  of  Sun  are  all  getting  closer  to 
those  of  SnU.  This  means  that  the  inclusion  size  effect  is  less  significant  and  may  be  ignored  for  large  inclusions  in  some 
cases,  which  agrees  with  the  general  trend  observed  experimentally  (e.g.,  Cho  et  al.,  2006). 

The  change  of  S1212  with  the  position  and  inclusion  size  is  illustrated  in  Figs.  9-12  together  with  a  comparison  with 
Sj212  for  the  three  types  of  polyhedral  inclusions.  Clearly,  S^2i2  varies  with  the  position  of  x  inside  each  polyhedral 
inclusion,  which  differs  from  that  in  an  ellipsoidal  inclusion  and  supports  the  Esheiby  conjecture.  But  the  classical  part 
Sj212  at  a  given  value  of  X]/R  remains  the  same  for  all  inclusion  sizes,  as  expected  from  the  discussion  in  Section  3.1. 
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Fig.  10.  Variation  of  S1212  along  the  x,  axis  inside  the  octahedral  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  edge  length  (see 
Fig.  5(b)). 


Xi/R  x-y/R  x1/R 

Fig.  11.  Variation  of  S1212  along  the  X\  axis  inside  the  tetrakaidecahedral  inclusion:  (a)  R=2L,  (b)  R=4L  and  (c)  R=6L,  with  R  being  half  of  the  cell  height 
(see  Fig.  5(c)). 


The  gradient  part  S^212,  as  the  difference  between  S1212  and  Sj212  (i.e.,  S^212  =  S12i2-S'f212),  is  seen  to  be  significantly  large 
for  small  inclusions  (e.g.,  R/L= 2)  and  becomes  insignificant  for  large  inclusions  for  all  three  types  of  polyhedral  inclusions. 
More  specifically,  it  is  observed  from  Fig.  9  that  for  the  cubic  inclusion  the  strain  gradient  effect,  as  measured  by  the  value 
°f  Sm2>  is  large  for  all  inclusion  sizes  when  x^/R>  0.6.  For  the  octahedral  inclusion,  the  strain  gradient  effect  is 
insignificant  and  can  be  neglected  when  R/L  >  4,  as  illustrated  in  Fig.  10.  For  the  tetrakaidecahedral  inclusion,  Fig.  1 1  shows 
that  the  strain  gradient  effect  is  also  small,  especially  in  the  region  away  from  the  square  faces. 

The  component  Sun  of  the  averaged  Eshelby  tensor  varying  with  the  inclusion  size  is  shown  in  Fig.  12  for  the  three 
inclusion  shapes,  where  Sllu  is  also  displayed  for  comparison.  The  values  of  Sun  shown  in  Fig.  12  are  obtained  using  Eqs. 
(41)  and  (39a),  (39b),  (39c)  and  (39d),  which  are  also  applied  to  get  the  values  of  S^nl  with  L-> 0. 

It  can  be  seen  from  Fig.  12  that  for  each  polyhedral  inclusion  Slnl  (based  on  classical  elasticity)  is  a  constant 
independent  of  the  inclusion  size  R.  However,  Snii  predicted  by  the  current  model  based  on  the  strain  gradient  elasticity 
theory  does  vary  with  the  inclusion  size:  the  smaller  the  inclusion,  the  smaller  the  Eshelby  tensor  component.  In 

particular,  when  the  inclusion  is  small,  the  strain  gradient  effect,  as  measured  by  Sml(  =  Snn-S^ln),  is  significantly  large 

—  — c 

and  should  not  be  ignored.  As  the  inclusion  becomes  large,  Smi  approaches  Sml  from  below,  indicating  that  the  strain 
gradient  effect  gets  small  and  may  be  neglected  for  very  large  inclusions. 

The  observations  made  here  are  also  true  for  the  other  components  of  the  Eshelby  tensor  S„w  in  Eqs.  (9a),  (37a),  (37b), 
(37c),  (37d),  (28)  and  (31)  and  its  volume  average  S,jW  in  Eq.  (41). 

From  the  numerical  results  presented  above,  it  is  clear  that  the  newly  obtained  Eshelby  tensor  based  on  the  SSGET  can 
capture  the  inclusion  size  effect  at  the  micron  scale,  while  the  Eshelby  tensor  based  on  classical  elasticity  does  not  have 
this  capability. 

5.  Conclusions 

An  analytical  solution  is  provided  for  the  Eshelby-type  problem  of  an  arbitrarily  shaped  polyhedral  inclusion  embedded 
in  an  infinite  elastic  matrix  using  a  simplified  strain  gradient  elasticity  theory  (SSGET)  that  contains  one  material  length 
scale  parameter  in  addition  to  two  classical  elastic  constants.  The  SSGET-based  Eshelby  tensor  for  a  polyhedral  inclusion  of 
arbitrary  shape  is  analytically  derived  in  a  general  form  in  terms  of  three  potential  functions,  two  of  which  are  the  same  as 
the  ones  involved  in  the  Eshelby  tensor  based  on  classical  elasticity.  These  potential  functions,  as  three  volume  integrals 
over  the  inclusion,  are  evaluated  by  dividing  the  polyhedral  inclusion  domain  into  tetrahedral  duplexes.  Each  of  the  three 
volume  integrals  is  first  transformed  to  a  surface  integral  by  applying  the  divergence  theorem,  which  is  then  transformed 
to  a  contour  (line)  integral  based  on  Stokes’  theorem  and  using  an  inverse  approach  that  differs  from  those  employed  in 
the  existing  studies  based  on  classical  elasticity.  The  newly  obtained  Eshelby  tensor  is  separated  into  a  classical  part  and  a 
gradient  part.  The  classical  part  depends  only  on  Poisson’s  ratio  of  the  matrix  material,  while  the  gradient  part  depends  on 
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Fig.  12.  Variation  of  Sun  with  the  inclusion  size:  (a)  cubic,  (b)  octahedral  and  (c)  tetrakaidecahedral. 


both  Poisson’s  ratio  and  the  material  length  scale  parameter  that  enables  the  explanation  of  the  inclusion  size  effect.  This 
SSGET-based  Eshelby  tensor  reduces  to  its  counterpart  based  on  classical  elasticity  when  the  strain  gradient  effect  is  not 
considered.  A  general  form  of  the  volume  averaged  Eshelby  tensor  over  the  polyhedral  inclusion  is  also  obtained,  which 
can  be  used  in  homogenization  analyses  of  composites  containing  polyhedral  inclusions. 

To  demonstrate  the  newly  derived  Eshelby  tensor,  three  types  of  polyhedral  inclusions,  cubic,  octahedral  and  tetra¬ 
kaidecahedral,  are  analyzed  by  directly  applying  the  general  formulas.  The  numerical  results  reveal  that  for  each  of  the  three 
inclusion  shapes  the  components  of  the  new  Eshelby  tensor  change  with  the  position  and  inclusion  size,  whereas  their  classical 
elasticity-based  counterparts  only  vary  with  the  position.  When  the  inclusion  is  small,  the  gradient  part  is  seen  to  contribute 
significantly  and  should  not  be  ignored.  Also,  it  is  found  that  the  smaller  the  inclusion  size  is,  the  smaller  the  components  of  the 
volume-averaged  Eshelby  tensor  are.  These  components  approach  from  below  the  values  of  their  classical  elasticity-based 
counterparts  as  the  inclusion  size  becomes  large.  Hence,  the  inclusion  size  effect  may  be  neglected  for  large  polyhedral 
inclusions  in  some  cases. 
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